!===========Programa que calcula determinante de matriz NxN============
!Programa requer um arquivo matriz_in.dat com a matriz, e NÃO pode
!existir o arquivo perm_out.dat

PROGRAM exer04

!======================================================================

IMPLICIT NONE

!=======================Declaração de variáveis========================

      INTEGER :: N, I, J, fat
      INTEGER, ALLOCATABLE, DIMENSION(:) :: vet
      REAL(4), ALLOCATABLE, DIMENSION(:,:) :: mat_in
      REAL(4) :: det, termo

!===================Cria os arquivos do modo desejado==================

      OPEN(UNIT=10, FILE='perm_in.dat', STATUS='unknown')

      !Escreve '1 1' no arquivo para poder gerar as permutações
      WRITE(10,*) '1 1'

      CLOSE(UNIT=10)

!==========Abertura de arquivos, leitura de dados e alocação===========

      OPEN(UNIT=10, FILE='matriz_in.dat', STATUS='old')

      WRITE(*,*) 'Entre com N (da matriz NxN)'
      READ(*,*) N

      ALLOCATE(mat_in(N,N), vet(N+1))

      READ(10,*) mat_in

!====================Calcula as permutações para N=====================

      CALL make_perm_out(N)

!==================Cálcula o determinante pela série===================

      !Abre o arquivo onde se encontram as permutações
      OPEN(UNIT=40, FILE='perm_out.dat', STATUS='old')

      !Inicializa det
      det = 0.0

      !Calculo do fatorial que é usado como limite no proximo loop
      fat = 1
      DO I=1,N
            fat = fat*I
      END DO

      !Vai lendo linha por linha do arquivo as permutações P_i e sua
      !paridade para com a ajuda delas, calcular o determinante por
      !uma série
      DO I=1,fat
            termo = 1.0
            !Lê a permutação P_i do arquivo
            READ(40,*) vet
            DO J=1,N
                  !Calcula o termo a ser somado na série
                  termo = termo*mat_in(J,vet(J))
            END DO
            !Soma o termo na série (sinal depende da paridade)
            det = det + vet(N+1)*termo
      END DO

      WRITE(*,*) det

!=========================Desaloca a memória===========================

      DEALLOCATE(mat_in,vet)

!======================Remove os arquivos gerados======================

      CALL SYSTEM('rm -f perm_out.dat')

!======================================================================

END PROGRAM exer04

!======================================================================
 


!=============================Subrotinas===============================

SUBROUTINE make_perm_out (N)

      IMPLICIT NONE

      INTEGER :: N, I, J, fat
      INTEGER, ALLOCATABLE, DIMENSION(:) :: leitura

      DO  I=1,(N-1)
            OPEN(UNIT=30, FILE='perm_in.dat', STATUS='old')
            OPEN(UNIT=20, FILE='perm_out.dat', STATUS='new', POSITION='append')

            ALLOCATE(leitura(I+1))

            !Calcula o fatorial que será usado como limite no proximo loop
            fat = 1
            DO J=1,I
                  fat = fat*J
            END DO
            !Lê linha por linha do arquivo de entrada e permuta a partir de
            !cada linha. Chama a subrotina para cada linha
            DO J=1,fat
                  READ(30,*) leitura
                  CALL permut(leitura,I)
            END DO
            
            CLOSE(UNIT=30)
            CLOSE(UNIT=20)

            DEALLOCATE(leitura)

            CALL SYSTEM ('mv perm_out.dat perm_in.dat')
      END DO

      !Muda o nome do arquivo para ficar como saída
      CALL SYSTEM ('mv perm_in.dat perm_out.dat')

END SUBROUTINE make_perm_out


!=================Subrotina que permuta dado o vetor===================

!Subrotina recebe o vetor no qual incluirá o próximo número e fará
!as permutações. Também recebe N (número inicial de valores fora a
!paridade)
SUBROUTINE permut (vet,N)

      IMPLICIT NONE

      INTEGER :: N, I, J, flag
      INTEGER, DIMENSION(N+1) :: vet
      INTEGER, DIMENSION(N+2) :: temp

      DO J=1,(N+1)
            !Inicializa o vetor temporário
            temp = 0
            !Verifica a paridade
            IF (vet(N+1) == 1) THEN
                  temp(N+2) = (-1)**(J+1)
            ELSE
                  temp(N+2) = (-1)**J
            END IF
            !Coloca o valor de N+1 na posição que depende de J
            temp(N+2-J) = (N+1)
            flag = 0

            !Loop que faz as permutações
            DO I=1,N
                  IF (temp(I) /= (N+1) .AND. flag == 0) THEN
                        !Para os valores diferentes e anteriores a N+1
                        temp(I) = vet(I)
                  ELSE
                        !Para os valores após o N+1
                        temp(I+1) = vet(I)
                        flag = 1
                  END IF
            END DO

            !Escreve no arquivo de saída as permutações geradas nesse loop
            WRITE(20,*) temp
      END DO

END SUBROUTINE permut

!======================================================================
